A new method to find full complex roots of a complex dispersion equation for light propagation 
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A new numerical method is presented to find full complex roots of a complex dispersion equation. For the 
application of the solution, the complex dispersion equation of a cylindrical metallic nanowire is investigated. 
By using this method, locus of Brewster angle, complex dispersion curves of Surface Plasmon Polaritons (SPPs) 
and complex bulk modes can be obtained in once calculation. Approximate analytical solution to the complex 
dispersion equation has also been derived to verify our method. 

PACS numbers: 73.20.Mf, 41.20.-q, 41.20.Jb, 78.20.Bh 



I. INTRODUCTION 

Dispersion relation is the basic property of light propa- 
gating in mediums, which specifies the relation between the 
wavevector k and the frequency a> of the light. Recently, 
the light propagation in metallic nanostructures has attracted 
enormous attentions from researchers due to the excitation of 
the Surface Plasmon-polaritons(SPPs) '"^ by the light. The 
SPPs are known for the wide application in nano optics at- 
tributed to the spatial localization of the SPPs at the metal- 
dielectric interfaces, which can be guided to manipulate light 
in nanoscaled photonic circuitry'"''. Furthermore, the spatial 
confinement of the SPPs can enhance the field intensity at 
the interfaces, which can influence the luminescence inten- 
sities^"** and life time of emitters close to the interfaces.'SiiS In 
order to understand the physical properties of SPPs in plas- 
monic structures, it is important to get the dispersion rela- 
tions of light propagation in the plasmonic structures. Princi- 
pally, the dispersion relations can be obtained by solving the 
Maxwell equations with the boundary conditions of the plas- 
monic structures imposed. For simple plasmonic structures, 
such as planar metal surfaces, the dispersion of SPPs can be 
solved analytically.""'^ But for the most structures which are 
nonsymmetrical and irregular, codes such as finite-element 
method (FEM) and finite-difference time-domain techniques 
(FDTD) need to be applied for numerical calculations. The 
basic technique of FEM or FDTD is to discrete the Maxwell 
equations, which are known for the time consumption. For- 
tunately, for some symmetric structures, such as cylindrical 



, ^ nanowires 



11,15-19 



the Maxwell equations imposed with the 



boundary conditions can be transformed into a dispersion 
equation. Solving the dispersion equation can improve the 
efficiency to get the dispersion relation rather than the numer- 
ical calculation by the discretion of Maxwell equations. 

The dispersion equation normally is nonlinear and tran- 
scendental. What is more, when the metal loss is introduced 
into the metal dielectric, '-^ '''^■^"'^' the dispersion equation then 
is a complex and transcendental equation. Thus, the disper- 
sion relations obtained from the complex equation are also 
complex. There exist two types of the complex dispersion re- 
lations for one same plasmonic structure, specifying either a 
complex frequency as a function of a real wave vector, noted 
as complex - co for convenience, or a complex wave vec- 
tor as a function of a real frequency, noted as complex - k. 
The two types of dispersion relations are rather different, even 



though they are for one same structure. The back bending as 
a characteristic of the complex - k relation is absent in the 
complex - (jj relation. The latter is an asymptotic curve."*' - 
It has been suggested that the complex - k solution of the 
dispersion relation describes the SPPs mode decaying spa- 
tially while the complex - oj solution is for the SPPs decay- 
ing in time rather than in space The discrepancy of 
the two solutions has been considered to be originated from 
the metal Ohmic loss i'^-'^'^"'^' Such conclusion is originally 
drawn from the study of SPPs on planar metallic-dielectric 
surfaces. When a perfect metal without damping is consid- 
ered in this planar case, the two types of solutions of the SPPs 
dispersion relations overlap with a same asymptotic behavior. 
When the metal Ohmic loss is introduced into the dielectric 
response of metal, the two types of solutions are then differ- 
ent. The original of the discrepancy of the two solutions due to 
the metal loss has also been confirmed in cylindrical metallic 
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The main mission to get the dispersion relation of the light 
propagation in plasmonic structures then is to find the com- 
plex roots of the equation. The complex dispersion equation 
can be described by f(x, y,z) - with three real variables x, 
y and z- For example, in the complex - o) solution, x, y and z 
represent the real part of oj (Re[a>]), the imaginary part of ll) 
{Im[a)]) and real wavevector k respectively. In the complex-k 
solution, X, y and z then represent the real part of k (Re[k]), 
imaginary part of k (Im[k]) and real oj respectively. To solve 
this complex equation, normally one variable is given, say z, 
then the complex equation can be simplified to be a complex 
equation f(x, y) = Q with only two real variables. One com- 
monly used method to solve this equation is to choose all pos- 
sible values of x and y in their given ranges to check if they sat- 
isfy the equation or not. To achieve this purpose, a coordina- 
tion system with x and y axis is gridded with an enough small 
mesh size and then the values of x and y at each grid point are 
substituted into the equation. Considering the error allowance 
6, the criterion for this grid method is to find the roots {xo,yo) 
if they satisfy \f(xo,yo)\ < 6. To this grid method, there exit 
two main short comings. The first is that in order to get the 
full solutions to the equation, one needs to gird the coordi- 
nate system with an enough small mesh size. Or, the solutions 
may be lost. However, in the reality for the calculation, we 
find that the mesh size can reach the order of the magnitude 
of 10"'^ and even smaller, which increases the computation 
time. The second shortcoming is that the criterion can not 
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FIG. 1: Schematic of our method to get full complex roots of a com- 
plex dispersion equation. 5 , and 5^, are the mesh sizes along x and y 
axis respectively, (a, b) and (c, d) are assumed to be the roots of the 
equation, enclosed in mesh A and mesh B respectively. The former is 
used for the linear approximation, while the latter is for the parabola 
approximation. 



guarantee the (x(),3'o) to be the root even though they satisfy 
the criterion since 6 is not rigorously equal to zero. No matter 
how small the 6 is taken, the second shortcoming still remains 
since 5 9^ 0. To improve the computation efficiency, the alter- 
native method is the Newton-Raphson (NR) method, which 
can converge to the roots of the equation quickly. However, it 
is well known that the NR method may miss the roots if there 
exist multiroots of the equation. What is more, the NR method 
may not converge if the initial estimate is not close enough to 
the root, or may converge to wrong root. Especially for the 
complex transcendental equation, which is very common for 
the light propagation in plasmonic structures when the metal 
Ohmic loss is introduced, the NR method will lose its power 
to find the roots. In this paper, we propose one new method to 
find full complex roots of complex transcendental equations. 
As an example for the application, the dispersion equation of 
light propagating in a cyUndrical metallic nanowire is investi- 
gated. 



II. METHOD 

A. linear approximation 

The main mission of our method to solve the complex equa- 
tion f{x, y) = is to find criterions used to judge if the val- 
ues of the variables are the roots or not. Similarly to the grid 
method, the coordination system with x and y axis is gridded 
firstly. Suppose that a and b are the roots of the equation sat- 
isfying f(a, b) -0 and the point (a, b) is enclosed in the mesh 
A(fig. [1]), then the function at each point of the mesh corner 
can be expanded at the point (a, b) by the first order Taylor 
series as 

f(x,yh = f(a,b) + f:(a,b)(-Aai)+f;(a,b)(Abi), (la) 
f(x,y)2^f(a,b)+f^(a,b)(Aa2) + f;(a,b)(Abi), (lb) 



f(x,yh = f(a, b) + f.ia, b){-Aa,) + f[.{a, b)(-Ab2), (Ic) 



f(x,y)4 = f(a, b) + b)(Aa2) + f;(a, b)(-Ab2). (Id) 

Here, subscript j of /(ji;,y)j(j=i,2,3,4) represents the calcu- 
lated results of the function f(x, y) at the corner point j{j = 
1,2, 3,4) of mesh A, which has been labeled in the fig.[T] The 
complex f'^fy^ia, b) is the partial derivative of the function with 
respect to x(y) at the point of {a,b). Aa\, Aa2, Ab\, and Ab2 
are the absolute values of the coordination components of the 
four corner points of the mesh A with respect to the point (a, b) 
considered as the coordinate origin. Eg. ([Tall subtracted from 
eq.lfTbli gives 

f'r h\ /(^^yh- f(x,y)i f(x,y)2-f(x,y)i 
fMb) = — — = , (2a) 

Afli + Aa2 i X 

and eq.lfTal) subtracted from eq.(fTc]i gives 

ft h\ f(x,y)i-f(x,yh fix,y)i - f(x,yh 
fM^^ = ^ aa " V ■ ^^^^ 

Abi + Ab2 Sy 

The S X and S , are the mesh sizes along the x and y axis re- 
spectively, which are given as constants when gridding the co- 
ordination system. We substitute eq.© into eq.lfTal. obtaining 



Aa\ — 



Abi - 



L2L12 - L1L22 
J-'iiLn - L11L22 



L\L2\ - L2L11 



(3a) 



(3b) 



with 



Li = Re[f(x,yh], L2 = Im[f(x,y)x], L„ = Re[-f,{a,b)], 
Ln^Re[f;.(a,b)], L21 = Im[-f:(a,b)], L22 = Im[f;(a,b)]. 



Similarly, by substituting eq.© into eq.lfTbb. ( flcb and ( fTdl i. 

we can get two sets of the Aoi, Afl2, Abi, and Ab2, noted as 
Aa"^""'fi,, and A/?"*""'^ ' for convenience. Now we give the 

m(m-l,2)' m{m=l,2) ^ 

criterion to find the roots of the equation. If the a and b are 
the roots of the complex equation satisfying /(a, b) - 0, then 



^ »(»-l,2) ^ ri 

Aa \ L > 0, 

m(m- 1,2) ' 



Aa 



1 

m(m— 1,2) 



Aa 



m(m= 1,2)' 



Abi > 0, 



m(m- 1 ,2) 



m(m- 1,2)' 



(4a) 



must be held. Or, there exists no root in the grid A due to 
the f(a,b) in eq.([TJ. Considering the error allowance 6, 
criterion (l4a]i can be modified as 



. «(«=1,2) ^ ry 

Aa , L > 0, 

m(m= 1,2) ' 



Abi > 0, 



<5, 



(4b) 



< 6. 



By using this method, the correct mesh enclosing the root can 
be found. And then we can further gird the correct mesh and 
consider the grid point with the minimum value of [/(x,}?)! in 
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the mesh as the root of the equation, at least they are very 
close to the roots. We note that in our method the criterion 
|/(xo,yo)l < <5 is not fatal for the determination of the roots. 
And the function values [/(x,}?)! at the points considered as 
the roots may be larger than the error allowance 6 if one wants 
to save computation time. The criterion (|4]l can guarantee the 
points close to the exact roots if the points are enclosed in the 
correct meshes satisfying the criterion (|4]i. In our method, all 
the meshes are independent to each other. Thus, all possible 
roots can be found, overcoming the problem met by the NR 
method. What is more, this method focuses on the finding 
of correct meshes enclosing the roots instead of the finding 
of correct points matching the roots, which saves much more 
computation time. 



B. parabola approximation 

Generally, if the derivative f'^^y^{x,y) in eq.([T]i may equal to 
zero when the function f(x, y) has a parabola-like shape, we 
can expand the function with the Taylor series up to order two 

fix, y) = f(c, d) + /u(c, <i)Ac + /i,,(c, d)/\d 

+flxx{c, d)AC^ + f2yy(c, d)Ad^ + f2xy(c, d)AcAd. 

(5) 

Here, fix(y) is the coefficient for the first order Taylor expan- 
sion and the f2x(y)x(y) is the coefficient for the second order. 
If the c and d are the roots of the complex equation, then 
f(c, t/) = is held. For illustration, point (c, d) is enclosed in 
the mesh B, shown in the fig.[T] Each point of mesh B used for 
the calculation has been labeled by j{j = 1,2,3,4, 5, 6, 7, 8) in 
the fig. 1 . Then, we can get 

hxy{c,d) = — , (6a) 



/«(«=^'2'(c,^)and4"=^'2'(c,rf): 



- /iy(c, d) - Ifiyyic, d)Ad\ + fixvic, d)Aci 
^ 4f(x,y)6-3fix,yh-f(x,yh 
S, 



(7b) 



The eq.dTab and ( fTbl i are complex. Each equation can be sep- 
arated into two equations with respect to the real and imag- 
inary parts of the equation. Now we consider three situa- 
tions. The first situation is /i t (c, d) = 0. The second situ- 
ation is fiy(c,d) = 0. And the third is both fi^(c,d) - 
and /iy(c, d) = 0. For the first situation, we can substitute 
eq.© into the eq.jTali and get two sets of Aci and Adi, noted 
as Ac^^''"^'^' and Ad'^^''^^'^\ And for the second situation, by 
substituting eq.© into eq.jTbji, we also obtain two sets ofAci 
and Adi, noted as Ac^'^''^'^''*' and Ad'!^''^^'*\ Last, for the third 



situation, we then can get four sets of the roots Ac 



and Ad' 



,p(p= 1,2,3,4) 



,Mp= 1,2,3,4) 



from the eq.O. Now we give the criterion to 
find roots of the equation. If the c and d are the roots of the 
complex equation satisfying f(c,d) - , then for the eq.® 
the criterion 
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< s, 



<s. 



(8a) 



must be held. And for the first situation of f\x{c, d) - 0, addi- 
tional criterion is 





Acj -ACj 1 




Acj 1 




Ad\ -Ad\ 




Ad\ 



<5, 



<6. 



(8b) 



For the second situation of fiy(c, d) — 0, the additional crite- 
rion is 







■ 


l\c\ -Ac\ 


/2yv(c, d), noted as 
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Ad^-Adf 
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<5, 



< 5. 



(8c) 



2(/(x,,v)i+/(A:,y);)-4/(x,y)5 



Si 



fl /„ _ 2(f(x,yh+f(x,y)i)-4f{x,y)e 



(6b) 



And for the third situation of both fix(c, d) - and /iy(c, d) - 
0, the additional criterions of eq.dSbji and dScl ) both should be 
held. 



and 



Si 



y2 _ 2(f(x,y)2 +f(x,yW-Af(x,yh 



V ■'2w 



Sf. 



(6c) 



Here, the subscript j of f{x, y)j specifies the function f(x, y) 
calculated at the point j of the mesh B. After some algebra, 
we get the following linear equation system: 



/u.(c, d) - 2f2xx(c, d)ACx + f2xy(c, d)Ad\ 

^ 4f(x,yh-3f(x,yh-f(x,y)2 

Sx 



(7a) 



III. CALCULATION AND DISCUSSION 

A. dispersion equation 

In order to verify our method, the complex dispersion equa- 
tion of SPPs on the planar metallic surface has been inves- 
tigated. The rigorous analytical solutions to the dispersion 
equationii^ in that case have confirmed the validity of our 
method. However, in this paper we focus our attention on the 
solutions to the dispersion equation of a cylindrical metallic 
nanowire. The nanowire has the shape with the radius of r 
and an infinite length in a medium. The electromagnetic field 
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can be expanded with cylindrical harmonics. By solving the 
Maxwell equation with the boundary conditions imposed, the 
dispersion equation can be obtained as the following transcen- 
dental equation: 

HlP(k^r)k^^ -J„iKtr)kl^ 

Hi'\k^r)k,k% -J„(k,ir)kokj^ ^ ^ 

k^rHi^^'{k^r)kok, -nk,Hl-\k^r)k, -k,irj'jk,ir)kah nkJ„(k,,r)ko 
-nk,Hi^\k^r) k^rHi^^' (k^r)ko nk,J„(Kir) -k,irJ'„(Kir)ki 

(9) 

which repeats the reported resultsJ^ii^ Here, H^^^ is the first 
kind of Hankel function with the order of integer n and /„ 
is the Bessel function of «th order The Hankel and Bessel 
functions with denote represent the first order differentiation. 
kj{j=o,\) is the wave vector with the value of kj = -^fejco/c, 
where Eyis the dielectric function and the subscript j labels 
the quantities outside the nanowire (/' = 0) or inside it (j = 1). 
c is the speed of light in vacuum. The dielectric function of 
the metal can be expressed as 

" I (10) 
a>{(i> + IT) 

where cjp is the buUc-plasmon frequency and r is the bulk elec- 
tron relaxation rate^^ which reflects the metal Ohmic loss. 
Eq.® then is a complex equation with t introduced into the 
dielectric. Eoc in eq.(fTO]i is a constant for the general descrip- 
tion of the dielectric function of the metal, k, is the component 
of the SPPs wave vector along the cylinder axial and the radial 

components of the wave vectors are defined as kyj - -^k^j - ki. 

We note that the value of kro should be chosen to guarantee the 
imaginary part of kro to be positive since the light intensity 
should be decaying away from the metal cylinder. The disper- 
sion relation between k- and u then can be obtained from the 
eq.(|9]l numerically. For the calculation, we renormalize the 
(jj and T by u)p. Wave vector components of kj^^cu) , fe, and 
ki j{j=o,\) are renormalized by ojp/c and r is renormalized by 
c/cDp. 

B. complex - to solution 

Fig.|2]shows the complex-io dispersion relations calculated 
from the eq.® by using our method. For the comparison be- 
tween the published results'^ and our results, the dielectrics 
take the values of eo = 5.3 and e^o - 9.6 with r - 0.005647 in 
the calculation. The renormalized radius r takes the value of 
r = 0.5 for the nanowire, which is corresponding to the real 
radius of 26.9nm. Four orders (n = 0, 1, 2, 3) of the dispersion 
relations are calculated. The left column of the figures is for 
the relation of Re[aj] and k^ and the right column is for the 
relation of /m[w] and In the calculation, the mesh size is 
0.01 and the error allowance 6 is set to be (5 = 0.1. It is shown 
that there exist two dispersion branches in the figures of « = 
and n = 1, shown in figs.|2ja)-(d). One branch is an asymp- 
totic curve with the frequency approaching the surface plas- 
mon frequency ojsp - Vfco/(eo + foo) ~ 0.8, which represents 
the SPPs surface wave. Another branch having the frequency 



above 1 is not a surface wave but identified as the locus of 
the Brewster angle. In our calculation, the SPPs dispersion 
curve are found by the using of the criterion (HI while the locus 
of the Brewster angle needs the criterion ([8]), reflecting that the 
relation between the k, and Re[a)] behaves parabola-like when 
it is close to the locus of the Brewster angle. For the order 
n > 1, the relation between the a> and k, is nearly dispersion- 
less and no locus of the Brewster angle can be observed any- 
more. The imaginary parts of the frequencies are very small, 
which have been shown in the right figures. Our calculated re- 
sults are coincident to the published results,'** confirming the 
validity of our method. Our results exhibit that two dispersion 
branches can be obtained in once calculation by our method, 
however, which can not be achieved by the NP method if the 
frequency range includes the two curves. 



C. complex — k solution 

We have repeated the reported results'' by using our 
method with eo = 1, eoo = 1, t = 0.039062 and n - \. How- 
ever, we only show our result with n - in fig. [3] since in 
this case an approximate analytical solution can be derived to 
verify our calculated results. The renormalized radius of the 
nanowire for the calculation is r = 2, which is correspond- 
ing to the real radius of 105«mJ2 By using our method, full 
solutions to the complex dispersion equation can be obtained, 
exhibiting many modes in fig.[3ja) and (b). Those modes can 
be classified into three, which have been shown in fig. Oc)- 
(h). Each figure in the left side responses for the relation of 
(A) and Relk,}, and the corresponding figure at its right side is 
for the relation of u) and Imlk^}. Dispersion curves in fig.^c) 
and (d) exhibit the locus of the Brewster angle. '^ Fig. Oe) 
and (f) represent the SPPs mode with a back bending in the 
curve. The back bending is induced by the metal loss.^** For a 
nanowire fabricated with perfect metal, the SPPs mode then is 
an asymptotic curve to infinity when the frequency approaches 
the (jjsp - 1/ V2. The third class of the modes in fig. |3lg) 
and (h) have an infinite mode number and the Re{h] of those 
modes are small compared to Iin[k^]. Those modes have been 
defined as bulk modes (BMDs)," but never been reported 
in this case. Therefore, there exists one question that if the 
BMDs are wrong solutions in our calculation. To answer this 
question, we derive an approximate solution to the complex 
dispersion equation in the following to verify our method. 

Before the derivation, three main characteristics of the cal- 
culated BMDs should be clarified. The first characteristic is 
that the dispersion curves of the BMDs are all nearly parallel 
to the a> axis, meaning that /m[fc,] of the BMDs are indepen- 
dent to ol>. The second is that the dispersion curves of the 
BMDs have a period of A Im[k^] xr ^ ji. Here, A Im[k^] is the 
interval between the curves. The third characteristic is that 
the lower dispersion curves with oj < co^p shift their phases 
by A Im[k-] x r x n/l with respect to the upper curves of 

CO > iOsp. 

For the derivation, it is reasonable for us to consider the 
limit case of Re[k,] — > and Im[k^] — » 00 based on fig. [3fg) 
and (h). Thus, the Bessel and Hankel functions in eq.(|9|l can 
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FIG. 2: complex-u) dispersion relations calculated by our method for a cylindrical metallic nanowire with eo = 5.3, eoo = 9.6 and t = 0.005647. 
The renormalized radius r = 0.5 is used for the calculation, which is corresponding to the real radius of 26.9nm of the nanowire. 
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FIG. 3: complex - k dispersion relations calculated by our method for a cylindrical metallic nanowire with 6o = \, = 1 and r = 0.039062. 
In the calculation, r = 2 and « = are used. The renormalized radius r = 2 is corresponding to the real radius of 105;!m of the nanowire. 



panded by k, as 



2kl - k] 
2ik, 



For simplicity, the metal is assumed to be perfect with t - 0. 
Substituting eqs.lfTTTiinto eq.®, we can get a simple equation 
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FIG. 4: complex - k dispersion relations calculated by our method for cylindrical metallic nanowires with various radius r. In the calculation, 
^ = 1, 600 = 1, « = and T = 0.039062 are used. 
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FIG. 5: complex - k dispersion relations calculated by our method for a cylindrical metallic nanowire with various orders n. In the calculation, 
= i, = i, r = 2 and r = 0.039062 are used. 



side of eq. (fT2l l. the complex k, can be solved analytically as 

(13) 



In eq.(fT3]). m is an integer. Due to the approximation in the 
derivation, the analytical solutions are not coincident to the 
calculated results of the BMDs. However, we find that in the 



Im{k,] 



KjA+lmn 
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analytical solutions the Im[k~] is independent to the w, which 
is the first characteristic of the BMDs we calculated. Sec- 
ondly, the analytical solution of Im[k,] shows that the solu- 
tions have a period of A Im[k,] x r - 2n. Considering the 
minus before y in the right side of eq.(fT2]l. the period is then 
equal to A Im[k,] x r = n , which is the second characteris- 
tic of the calculated BMDs results. Last, for the upper curves 
with w > 1 / V2, the denominator of 7 is a pure imaginary 
value, which is equivalent to the phase shift of A Im[k,] x r 
by njl with respect to the lower curves. This conclusion is 
just the third characteristic of the BMDs results calculated by 
our method. Thus, the approximate analytical solution to the 
eq.(|9]) confirms the validity of our method. As shown in fig. 3, 
full solutions to the complex dispersion equation can be ob- 
tained in once calculation by our proposed method. 

In the derivation, r = means that the BMDs are induced 
by the structure formation instead of metal loss. For details, 
fig. E] shows the dispersion curves of BMDs with various r 
in the case of n = 0, eo = 1, ec« = 1, and t - 0.039062. 
The curves have two peaks with peak one close to the surface 
plamson frequency and peak two at low frequency. For the 
nanowires with a smaller radius r, the A Imlk-} is larger to 
hold the period A Im{kJ\ x r - n. In this case, peak one is 



larger while peak two smaller We also find that the period 
is the basic property of the BMDs, which is independent to 
the integer order n. The BMDs for different orders has been 
presented in fig.|5] showing the Im{k~] of BMDs have the same 
period but different position. We suggest that the order n only 
influences the initial position of Im{kJ\. 



IV. SUMMARY 

We have proposed one new method to find full complex 
roots of a complex transcendental equation. The correct 
meshes enclosing the roots are independent to each other, 
which guarantee the finding of the all roots. For the appli- 
cation of this method, the complex dispersion equation of a 
cylindrical metallic nanowire is investigated. In our calcula- 
tion, locus of the Brewster angle, SPPs dispersion curves, and 
bulk modes all can be obtained in once calculation. Approx- 
imate analytical solution to the dispersion equation has been 
derived to verify our results. This method can be applied to 
all other complex transcendental equations with two real vari- 
ables. 
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